Load and install R packages:
packages <- c('lmerTest', 'lme4', 'ggplot2', 'tidyverse', 'readxl', 'purrr', 'performance', 'emmeans', 'MASS', 'dplyr', 'tidyr')
# Check if each package is installed; if not, install it
for (pkg in packages) {
if (!requireNamespace(pkg, quietly = TRUE)) {
install.packages(pkg)
}
library(pkg, character.only = TRUE)
}
## Loading required package: lme4
## Loading required package: Matrix
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
##
##
## Attaching package: 'MASS'
##
##
## The following object is masked from 'package:dplyr':
##
## select
Grab neophobia_data.csv from processed_data
directory:
script_dir <- dirname(rstudioapi::getActiveDocumentContext()$path)
# Set the working directory to the script's directory
setwd(script_dir)
data <- read.csv("processed_data/neophobia_data.csv", row.names = 1)
Summary statistics of the data:
# Quick look at structure
str(data)
## 'data.frame': 536 obs. of 15 variables:
## $ Bird_ID : chr "BB_BB" "BB_BB" "BB_BB" "BB_BB" ...
## $ Trial : int 1 2 3 4 5 6 7 8 1 2 ...
## $ Enclosure : chr "B8" "B8" "B8" "B8" ...
## $ Object : chr "novel" "novel" "control" "control" ...
## $ Context : chr "group" "individual" "group" "individual" ...
## $ Latency_to_enter: num 2.2 6.23 2.27 1.57 1.07 ...
## $ Latency_to_Eat : num 2.9 13.83 2.2 3.73 2.73 ...
## $ Zoi_duration : num 183.7 91.9 85.4 203.3 371.5 ...
## $ Object_Type : int 2 4 3 3 1 5 3 3 5 1 ...
## $ GroupID : chr "7A" NA "7A" NA ...
## $ NestID : chr "2024VM08" "2024VM08" "2024VM08" "2024VM08" ...
## $ group_dummy : int 1 0 1 0 1 0 1 0 1 0 ...
## $ ind_dummy : int 0 1 0 1 0 1 0 1 0 1 ...
## $ Object_contrast : num 0.5 0.5 -0.5 -0.5 0.5 0.5 -0.5 -0.5 0.5 0.5 ...
## $ Context_contrast: num 0.5 -0.5 0.5 -0.5 0.5 -0.5 0.5 -0.5 0.5 -0.5 ...
summary(data)
## Bird_ID Trial Enclosure Object
## Length:536 Min. :1.00 Length:536 Length:536
## Class :character 1st Qu.:2.75 Class :character Class :character
## Mode :character Median :4.50 Mode :character Mode :character
## Mean :4.50
## 3rd Qu.:6.25
## Max. :8.00
## Context Latency_to_enter Latency_to_Eat Zoi_duration
## Length:536 Min. : 0.800 Min. : 1.792 Min. : 0.00
## Class :character 1st Qu.: 1.467 1st Qu.: 2.800 1st Qu.: 38.93
## Mode :character Median : 1.880 Median : 3.633 Median : 76.35
## Mean : 8.380 Mean : 42.457 Mean :128.20
## 3rd Qu.: 2.600 3rd Qu.: 5.684 3rd Qu.:162.03
## Max. :600.000 Max. :600.000 Max. :598.13
## Object_Type GroupID NestID group_dummy
## Min. :1.000 Length:536 Length:536 Min. :0.0
## 1st Qu.:2.000 Class :character Class :character 1st Qu.:0.0
## Median :3.000 Mode :character Mode :character Median :0.5
## Mean :2.787 Mean :0.5
## 3rd Qu.:4.000 3rd Qu.:1.0
## Max. :5.000 Max. :1.0
## ind_dummy Object_contrast Context_contrast
## Min. :0.0 Min. :-0.5 Min. :-0.5
## 1st Qu.:0.0 1st Qu.:-0.5 1st Qu.:-0.5
## Median :0.5 Median : 0.0 Median : 0.0
## Mean :0.5 Mean : 0.0 Mean : 0.0
## 3rd Qu.:1.0 3rd Qu.: 0.5 3rd Qu.: 0.5
## Max. :1.0 Max. : 0.5 Max. : 0.5
The GroupID values are currently coded as
NA during individual trials. Assign the most frequent group
for each bird:
most_frequent_group <- function(group_ids) {
group_ids <- group_ids[!is.na(group_ids)]
if (length(group_ids) == 0) return(NA)
names(sort(table(group_ids), decreasing = TRUE))[1]
}
most_frequent_group_per_bird <- data %>%
group_by(Bird_ID) %>%
summarize(most_frequent_group = most_frequent_group(GroupID)) %>%
ungroup()
data <- data %>%
left_join(most_frequent_group_per_bird, by = "Bird_ID") %>%
mutate(GroupID = ifelse(is.na(GroupID), most_frequent_group, GroupID)) %>%
dplyr::select(-most_frequent_group)
Adjust the Trial variable by subtracting 1:
data$Trial <- data$Trial - 1
The full model as described in the RR, includes
Object_contrast, Context_contrast, and
Trial, and a complex random effect structure:
enter_model <- lmer(Latency_to_enter ~ Object_contrast * Context_contrast + Trial +
(1 | NestID) +
(-1 + group_dummy | GroupID) +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data)
## boundary (singular) fit: see help('isSingular')
summary(enter_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Latency_to_enter ~ Object_contrast * Context_contrast + Trial +
## (1 | NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy +
## group_dummy | Bird_ID)
## Data: data
##
## REML criterion at convergence: 5712.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4529 -0.1408 -0.0422 0.0686 10.2384
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Bird_ID ind_dummy 1.198e+03 34.6124
## group_dummy 2.852e-02 0.1689 1.00
## NestID (Intercept) 0.000e+00 0.0000
## GroupID group_dummy 0.000e+00 0.0000
## Residual 2.268e+03 47.6258
## Number of obs: 536, groups: Bird_ID, 67; NestID, 43; GroupID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 17.4867 4.3340 279.0124 4.035 7.06e-05
## Object_contrast 5.1073 4.1143 464.9955 1.241 0.21510
## Context_contrast -12.0511 5.8852 79.6165 -2.048 0.04389
## Trial -2.6020 0.9052 471.1684 -2.874 0.00423
## Object_contrast:Context_contrast -9.7670 8.2292 464.9964 -1.187 0.23588
##
## (Intercept) ***
## Object_contrast
## Context_contrast *
## Trial **
## Object_contrast:Context_contrast
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Objct_ Cntxt_ Trial
## Objct_cntrs 0.005
## Cntxt_cntrs -0.356 0.000
## Trial -0.731 -0.007 0.007
## Objct_cn:C_ -0.010 0.000 0.000 0.013
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
check_model(enter_model)
The random effect structure seems to be too complex for the amount of
data. However, the variance for the (Intercept) under
NestID is 0 and the variance for
group_dummy under GroupID is 0.
This suggests that both nest as differences between groups contribute
minimally to the variance in the outcome model. Let’s drop both
effects.
enter_model2 <- lmer(Latency_to_enter ~ Object_contrast * Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data)
## boundary (singular) fit: see help('isSingular')
summary(enter_model2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Latency_to_enter ~ Object_contrast * Context_contrast + Trial +
## (-1 + ind_dummy + group_dummy | Bird_ID)
## Data: data
##
## REML criterion at convergence: 5712.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4529 -0.1408 -0.0422 0.0686 10.2384
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Bird_ID ind_dummy 1198.0356 34.6127
## group_dummy 0.0286 0.1691 1.00
## Residual 2268.2114 47.6257
## Number of obs: 536, groups: Bird_ID, 67
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 17.4867 4.3340 279.0054 4.035 7.06e-05
## Object_contrast 5.1073 4.1143 464.9960 1.241 0.21510
## Context_contrast -12.0511 5.8852 79.6162 -2.048 0.04389
## Trial -2.6020 0.9052 471.1689 -2.874 0.00423
## Object_contrast:Context_contrast -9.7670 8.2292 464.9969 -1.187 0.23588
##
## (Intercept) ***
## Object_contrast
## Context_contrast *
## Trial **
## Object_contrast:Context_contrast
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Objct_ Cntxt_ Trial
## Objct_cntrs 0.005
## Cntxt_cntrs -0.356 0.000
## Trial -0.731 -0.007 0.007
## Objct_cn:C_ -0.010 0.000 0.000 0.013
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
check_model(enter_model2)
Seems like the interaction between object and context is non-significant, let’s drop it. Given the non-normal distribution, let’s logtransform the data as well:
enter_model3 <- lmer(log(Latency_to_enter) ~ Object_contrast + Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data)
## boundary (singular) fit: see help('isSingular')
summary(enter_model3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Latency_to_enter) ~ Object_contrast + Context_contrast +
## Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
## Data: data
##
## REML criterion at convergence: 1123.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9335 -0.4596 -0.0811 0.3538 6.9301
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Bird_ID ind_dummy 0.307211 0.55427
## group_dummy 0.008857 0.09411 1.00
## Residual 0.386408 0.62162
## Number of obs: 536, groups: Bird_ID, 67
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.29304 0.06319 188.76075 20.464 <2e-16 ***
## Object_contrast 0.10177 0.05370 465.99922 1.895 0.0587 .
## Context_contrast -0.19502 0.07774 76.68931 -2.508 0.0142 *
## Trial -0.13771 0.01179 468.78462 -11.679 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Objct_ Cntxt_
## Objct_cntrs 0.004
## Cntxt_cntrs -0.458 0.000
## Trial -0.653 -0.007 0.007
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
check_model(enter_model3)
There are still issues with the data distribution, let’s remove birds that did not interact and boxcox_transform:
data_enter <- data %>% filter(Latency_to_enter != 600)
boxcox_transform <- boxcox(lm(Latency_to_enter ~ Object_contrast + Context_contrast + Trial, data = data_enter))
best_lambda <- boxcox_transform$x[which.max(boxcox_transform$y)]
data_enter$Latency_to_enter_trans <- (data_enter$Latency_to_enter^best_lambda - 1) / best_lambda
enter_model4_boxcox <- lmer(Latency_to_enter_trans ~ Object_contrast +
Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_enter)
summary(enter_model4_boxcox)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Latency_to_enter_trans ~ Object_contrast + Context_contrast +
## Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
## Data: data_enter
##
## REML criterion at convergence: -115.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.00860 -0.62234 0.05268 0.62220 2.42212
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Bird_ID ind_dummy 0.023257 0.15250
## group_dummy 0.007021 0.08379 0.84
## Residual 0.037006 0.19237
## Number of obs: 533, groups: Bird_ID, 67
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.686462 0.020720 165.959664 33.130 <2e-16 ***
## Object_contrast 0.007148 0.016674 395.856100 0.429 0.668
## Context_contrast -0.028455 0.020248 64.608541 -1.405 0.165
## Trial -0.051896 0.003669 424.210940 -14.146 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Objct_ Cntxt_
## Objct_cntrs 0.007
## Cntxt_cntrs -0.301 -0.002
## Trial -0.623 -0.009 0.016
check_model(enter_model4_boxcox)
Comparing the different models, the boxcox transfromed comes out as best:
enter_model <- lmer(Latency_to_enter ~ Object_contrast * Context_contrast + Trial +
(1 | NestID) +
(-1 + group_dummy | GroupID) +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_enter)
## boundary (singular) fit: see help('isSingular')
enter_model2 <- lmer(Latency_to_enter ~ Object_contrast * Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_enter)
## boundary (singular) fit: see help('isSingular')
enter_model3 <- lmer(log(Latency_to_enter) ~ Object_contrast + Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_enter)
## boundary (singular) fit: see help('isSingular')
enter_model4_boxcox <- lmer(Latency_to_enter_trans ~ Object_contrast +
Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_enter)
anova(enter_model, enter_model2, enter_model3, enter_model4_boxcox)
## refitting model(s) with ML (instead of REML)
## Data: data_enter
## Models:
## enter_model3: log(Latency_to_enter) ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
## enter_model4_boxcox: Latency_to_enter_trans ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
## enter_model2: Latency_to_enter ~ Object_contrast * Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
## enter_model: Latency_to_enter ~ Object_contrast * Context_contrast + Trial + (1 | NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy + group_dummy | Bird_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## enter_model3 8 959.8 994.0 -471.91 943.8
## enter_model4_boxcox 8 -127.5 -93.2 71.73 -143.5 1087.3 0
## enter_model2 9 5148.8 5187.3 -2565.41 5130.8 0.0 1 1
## enter_model 11 5152.8 5199.9 -2565.41 5130.8 0.0 2 1
A similar model is built for Latency_to_Eat, again
incorporating interaction terms and random effects:
eat_model <- lmer(log(Latency_to_Eat) ~ Object_contrast * Context_contrast + Trial +
(1 | NestID) +
(- 1 + group_dummy | GroupID) +
(- 1 + ind_dummy + group_dummy | Bird_ID),
data = data)
## boundary (singular) fit: see help('isSingular')
summary(eat_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Latency_to_Eat) ~ Object_contrast * Context_contrast + Trial +
## (1 | NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy +
## group_dummy | Bird_ID)
## Data: data
##
## REML criterion at convergence: 1636.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1780 -0.5149 -0.1108 0.1927 4.6862
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Bird_ID ind_dummy 0.78398 0.8854
## group_dummy 0.08989 0.2998 1.00
## NestID (Intercept) 0.00000 0.0000
## GroupID group_dummy 0.00000 0.0000
## Residual 1.00669 1.0033
## Number of obs: 536, groups: Bird_ID, 67; NestID, 43; GroupID, 15
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 1.89202 0.10740 165.73052 17.616
## Object_contrast 0.74392 0.08668 465.00491 8.583
## Context_contrast -1.00808 0.11239 90.13817 -8.969
## Trial -0.04274 0.01899 466.49134 -2.251
## Object_contrast:Context_contrast -0.85108 0.17337 465.00516 -4.909
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Object_contrast < 2e-16 ***
## Context_contrast 3.94e-14 ***
## Trial 0.0248 *
## Object_contrast:Context_contrast 1.27e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Objct_ Cntxt_ Trial
## Objct_cntrs 0.004
## Cntxt_cntrs -0.434 0.000
## Trial -0.619 -0.007 0.008
## Objct_cn:C_ -0.008 0.000 0.000 0.013
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
check_model(eat_model)
Similarly, the variance for the (Intercept) under
NestID is 0 and the variance for
group_dummy under GroupID is 0.
Let’s drop those:
eat_model2 <- lmer(log(Latency_to_Eat) ~ Object_contrast * Context_contrast + Trial +
(- 1 + ind_dummy + group_dummy | Bird_ID),
data = data)
## boundary (singular) fit: see help('isSingular')
summary(eat_model2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Latency_to_Eat) ~ Object_contrast * Context_contrast + Trial +
## (-1 + ind_dummy + group_dummy | Bird_ID)
## Data: data
##
## REML criterion at convergence: 1636.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1779 -0.5149 -0.1108 0.1927 4.6862
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Bird_ID ind_dummy 0.78392 0.8854
## group_dummy 0.08988 0.2998 1.00
## Residual 1.00669 1.0033
## Number of obs: 536, groups: Bird_ID, 67
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 1.89202 0.10740 165.74889 17.616
## Object_contrast 0.74392 0.08668 465.00248 8.583
## Context_contrast -1.00808 0.11239 90.14508 -8.969
## Trial -0.04274 0.01899 466.48897 -2.251
## Object_contrast:Context_contrast -0.85108 0.17337 465.00273 -4.909
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Object_contrast < 2e-16 ***
## Context_contrast 3.94e-14 ***
## Trial 0.0248 *
## Object_contrast:Context_contrast 1.27e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Objct_ Cntxt_ Trial
## Objct_cntrs 0.004
## Cntxt_cntrs -0.434 0.000
## Trial -0.619 -0.007 0.008
## Objct_cn:C_ -0.008 0.000 0.000 0.013
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
check_model(eat_model2)
There are still some issues with the normality. Let’s remove non participating birds, and try the same boxcox tranformation:
data_eat <- data %>% filter(Latency_to_Eat != 600)
boxcox_transform <- boxcox(lm(Latency_to_Eat ~ Object_contrast * Context_contrast + Trial, data = data))
best_lambda <- boxcox_transform$x[which.max(boxcox_transform$y)]
data_eat$Latency_to_eat_trans <- (data_eat$Latency_to_Eat^best_lambda - 1) / best_lambda
enter_model4_boxcox <- lmer(Latency_to_eat_trans ~ Object_contrast *
Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_eat)
summary(enter_model4_boxcox)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Latency_to_eat_trans ~ Object_contrast * Context_contrast + Trial +
## (-1 + ind_dummy + group_dummy | Bird_ID)
## Data: data_eat
##
## REML criterion at convergence: -772
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3298 -0.6270 -0.0665 0.5674 3.4392
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Bird_ID ind_dummy 0.007875 0.08874
## group_dummy 0.005565 0.07460 0.92
## Residual 0.009093 0.09536
## Number of obs: 506, groups: Bird_ID, 67
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 0.782688 0.012637 123.643952 61.936
## Object_contrast 0.058827 0.008573 375.871518 6.862
## Context_contrast -0.153939 0.009588 63.404669 -16.056
## Trial 0.001207 0.001879 399.950402 0.642
## Object_contrast:Context_contrast -0.013930 0.017115 377.642252 -0.814
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Object_contrast 2.81e-11 ***
## Context_contrast < 2e-16 ***
## Trial 0.521
## Object_contrast:Context_contrast 0.416
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Objct_ Cntxt_ Trial
## Objct_cntrs 0.042
## Cntxt_cntrs -0.186 -0.050
## Trial -0.534 -0.039 0.042
## Objct_cn:C_ -0.041 -0.067 0.056 0.041
check_model(enter_model4_boxcox)
# now without interaction as it's non-sign:
eat_model4_boxcox <- lmer(Latency_to_eat_trans ~ Object_contrast +
Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_eat)
summary(eat_model4_boxcox)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Latency_to_eat_trans ~ Object_contrast + Context_contrast + Trial +
## (-1 + ind_dummy + group_dummy | Bird_ID)
## Data: data_eat
##
## REML criterion at convergence: -777.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3645 -0.6266 -0.0550 0.5560 3.4109
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Bird_ID ind_dummy 0.007853 0.08862
## group_dummy 0.005569 0.07463 0.92
## Residual 0.009088 0.09533
## Number of obs: 506, groups: Bird_ID, 67
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.782267 0.012621 123.342323 61.981 < 2e-16 ***
## Object_contrast 0.058358 0.008551 376.056191 6.825 3.54e-11 ***
## Context_contrast -0.153497 0.009564 63.210493 -16.049 < 2e-16 ***
## Trial 0.001270 0.001877 401.021902 0.677 0.499
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Objct_ Cntxt_
## Objct_cntrs 0.040
## Cntxt_cntrs -0.182 -0.047
## Trial -0.533 -0.036 0.040
check_model(eat_model4_boxcox)
boxcox model seems to fit data best
eat_model <- lmer(Latency_to_Eat ~ Object_contrast * Context_contrast + Trial +
(1 | NestID) +
(-1 + group_dummy | GroupID) +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_eat)
## boundary (singular) fit: see help('isSingular')
eat_model2 <- lmer(Latency_to_Eat ~ Object_contrast * Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_eat)
## boundary (singular) fit: see help('isSingular')
eat_model3 <- lmer(log(Latency_to_Eat) ~ Object_contrast + Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_eat)
## boundary (singular) fit: see help('isSingular')
eat_model4_boxcox <- lmer(Latency_to_eat_trans ~ Object_contrast +
Context_contrast + Trial +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data_eat)
anova(eat_model, eat_model2, eat_model3, eat_model4_boxcox)
## refitting model(s) with ML (instead of REML)
## Data: data_eat
## Models:
## eat_model3: log(Latency_to_Eat) ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
## eat_model4_boxcox: Latency_to_eat_trans ~ Object_contrast + Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
## eat_model2: Latency_to_Eat ~ Object_contrast * Context_contrast + Trial + (-1 + ind_dummy + group_dummy | Bird_ID)
## eat_model: Latency_to_Eat ~ Object_contrast * Context_contrast + Trial + (1 | NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy + group_dummy | Bird_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## eat_model3 8 1024.9 1058.7 -504.45 1008.9
## eat_model4_boxcox 8 -794.8 -761.0 405.39 -810.8 1819.7 0
## eat_model2 9 5070.9 5108.9 -2526.44 5052.9 0.0 1 1
## eat_model 11 5074.9 5121.4 -2526.44 5052.9 0.0 2 1
The full model as described in the RR, includes
Object_contrast, Context_contrast, and
Trial, and a complex random effect structure:
zoi_model <- lmer(Zoi_duration ~ Object_contrast * Context_contrast + Trial +
(1 | NestID) +
(-1 + group_dummy | GroupID) +
(-1 + ind_dummy + group_dummy | Bird_ID),
data = data)
summary(zoi_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Zoi_duration ~ Object_contrast * Context_contrast + Trial + (1 |
## NestID) + (-1 + group_dummy | GroupID) + (-1 + ind_dummy +
## group_dummy | Bird_ID)
## Data: data
##
## REML criterion at convergence: 6621.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9343 -0.5729 -0.2174 0.2772 3.5911
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Bird_ID ind_dummy 5689 75.43
## group_dummy 1638 40.47 0.81
## NestID (Intercept) 868 29.46
## GroupID group_dummy 1495 38.66
## Residual 11668 108.02
## Number of obs: 536, groups: Bird_ID, 67; NestID, 43; GroupID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 113.131 12.909 89.748 8.764 1.09e-13
## Object_contrast 3.697 9.332 399.031 0.396 0.692
## Context_contrast 17.667 14.913 20.537 1.185 0.250
## Trial 4.116 2.066 412.595 1.992 0.047
## Object_contrast:Context_contrast 106.889 18.664 399.033 5.727 2.01e-08
##
## (Intercept) ***
## Object_contrast
## Context_contrast
## Trial *
## Object_contrast:Context_contrast ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Objct_ Cntxt_ Trial
## Objct_cntrs 0.004
## Cntxt_cntrs 0.099 0.000
## Trial -0.559 -0.007 0.007
## Objct_cn:C_ -0.007 0.000 0.000 0.013
check_model(zoi_model)